Analysis CLL Proteomics - Final data Limma Proteomics
Sophie Rabe
Analysis date: 2020-01-10
Setup
Load libraries
library(plyr)
library(gtools)
library(openxlsx)
library(pheatmap)
library(reshape2)
library(progress)
library(Matrix)
library(Hmisc)
library(lemon)
library(ggpubr)
library(effsize)
library(ggbeeswarm)
library(ggfortify)
library(ggpmisc)
library(ggrepel)
library(readxl)
library(DESeq2)
library(TOSTER)
library(tidyverse)
library(vsn)
library(fdrtool)
library(limma)
library(apeglm)
library(IHW)
library(Rtsne)
library(biomartr)
library(biomaRt)
library(MultiAssayExperiment)
library(PMA)
library(gplots)
library(RColorBrewer)
library(grid)
library(ConsensusClusterPlus)
library(survival)
library(survminer)
library(cowplot)
library(matrixStats)
library(DEqMS)Load data
source("/Volumes/sd17b003/Sophie/Analysis/Screen_analysis/Figure_layouts.R")
load("/Volumes/sd17b003/Sophie/Analysis/CLL_Proteomics/CLL_Proteomics_final/Proteomics_Git/Robjects/CLL_Proteomics_Setup.RData")
load("/Volumes/sd17b003/Sophie/Analysis/CLL_Proteomics/CLL_Proteomics_final/Proteomics_Git/Robjects/CLL_Proteomics_ConsensusClustering.RData")colData(multiomics_MAE)$PG <- as.factor(CCP_group5[rownames(colData(multiomics_MAE))])
colData(multiomics_MAE)$CCP6_RNA <- as.factor(CCP_group6_RNA[rownames(colData(multiomics_MAE))])Limma Proteomics and genetic alterations
Functions
Function with DEQMS
Calculate_Limma <- function(assay_data, mut){
# prepare assay data
assay_data[is.infinite((assay_data))] <- NA
is.na(assay_data) <- NA
stopifnot(mut %in% colnames(proteomics_tbl_meta_biomart))
# prepare row data
row_data <- left_join((rownames(assay_data) %>% enframe(value = "rowname")),
(proteomics_tbl_meta_biomart %>%
dplyr::select(rowname, start_position, end_position, chromosome_name, mean_position) %>% unique %>%
group_by(rowname) %>%
dplyr::slice(1) %>% ungroup() ),
by=c("rowname")) %>% dplyr::select(-name) %>% as.data.frame()
rownames(row_data) <- row_data$rowname
stopifnot(all(rownames(assay_data)==rownames(row_data) ))
# prepare column data
col_data <- left_join((colnames(assay_data) %>% enframe(value = "colname")),
(proteomics_tbl_meta_biomart %>% dplyr::select(colname, mut) %>% unique() ),
by=c("colname") ) %>% dplyr::select(-name) %>% as.data.frame()
rownames(col_data) <- col_data$colname
stopifnot(all(colnames(assay_data)==rownames(col_data) ))
if( (length(unique(col_data[,mut])) ==2) & !is.logical(col_data[,mut]) & is.numeric(col_data[,mut]) ){
col_data[,mut] <- as.logical(col_data[,mut] == max(col_data[,mut]) )
}
#Creating an expression set
raw_dataE <- ExpressionSet(assayData = assay_data,
phenoData = AnnotatedDataFrame(col_data),
featureData = AnnotatedDataFrame(row_data))
validObject(raw_dataE)
if(!is.logical(pData(raw_dataE)[,colnames(pData(raw_dataE))==mut]))stop("The mutation you want to look at is not logical")
raw_dataE
# Limma
limma_data <- raw_dataE
comparison <- c("mut - wt")
colnames(pData(limma_data))[colnames(pData(limma_data)) == mut] <- "condition"
pData(limma_data) <- mutate(pData(limma_data),condition= if_else(condition, "mut","wt") )
limma_data <- limma_data[, !is.na(pData(limma_data)$cond)]
limma.cond <- factor(pData(limma_data)$condition, ordered = FALSE)
contrast.matrix <- model.matrix( ~ 0 + limma.cond)
colnames(contrast.matrix) <- gsub("limma.cond", "", colnames(contrast.matrix))
limma.object <- eBayes(
contrasts.fit(
lmFit(limma_data, design = contrast.matrix),
makeContrasts(contrasts = comparison, levels = contrast.matrix)
)
)
##### DEQMS
tmp <- metadata(multiomics_MAE)$protein_description %>% dplyr::select(`Gene Name`, MinPSMQuant) %>%
mutate(MinPSMQuant= as.numeric(MinPSMQuant))
psm.count.table <- tmp[,2] %>% as.data.frame()
rownames(psm.count.table) <- tmp$`Gene Name`
limma.object$count = psm.count.table[rownames(limma.object$coefficients),"MinPSMQuant"]
limma_DEQMS = spectraCounteBayes(limma.object)
DEqMS.results = outputResult(limma_DEQMS,coef_col = 1)
limma_results_i <- DEqMS.results
limma_results_i <- subset(limma_results_i, !is.na(logFC))
limma_results_i$comparison <- comparison
limma_results_i$mut <- mut
colnames(limma_results_i)[c(8,9,10, 14, 15, 16)] <-
c("limma.t_beforeDEQMS", "pvalue.limma_beforeDEQMS" , "fdr.limma_beforeDEQMS", "t", "pvalue.limma", "fdr.limma")
limma_results_i$fdr <- limma_results_i$fdr.limma
limma_results_i$pvalue <- limma_results_i$pvalue.limma
return(limma_results_i)
}Annotation function
Annotate_Limma_Results <- function(limma_results){
fdr_hit_threshold <- 0.001
fdr_candidate_threshold = 0.01
fc_hit_threshold <- log2(1.5)
fc_candidate_threshold <- log2(1.2)
limma_results$hit <-
with(limma_results, ifelse(fdr <= fdr_hit_threshold & abs(logFC) >= log2(fc_hit_threshold), TRUE, FALSE))
limma_results$hit_annotation <- with(limma_results,
ifelse(fdr <= fdr_hit_threshold & abs(logFC) >= log2(fc_hit_threshold),
"hit",
ifelse(fdr <= fdr_candidate_threshold & abs(logFC) >= log2(fc_candidate_threshold),
"candidate", "no hit")))
limma_results$hit_annotation <- factor(limma_results$hit_annotation, ordered = TRUE, levels = c("hit", "candidate", "no hit"))
return(limma_results)
}All genetic alterations limma
Conduct limma
limma_results <- NULL
mut_to_limma <- colnames(proteomics_tbl_meta_biomart)[-(1:11)][c(-15)]
for(i in 1:length(mut_to_limma)){
m <- mut_to_limma[i]
print(m)
limma_results <- bind_rows(limma_results,
Calculate_Limma(assay_data = assay(multiomics_MAE[prot_few_nas ,pat_overlap_prot_RNA ,"proteomics"]), mut = m))
}## [1] "chrom_abber_del11q"
## [1] "chrom_abber_del13q14"
## [1] "chrom_abber_del17p13"
## [1] "chrom_abber_gain8q24"
## [1] "chrom_abber_trisomy12"
## [1] "SNPs_ATM"
## [1] "SNPs_BIRC3"
## [1] "SNPs_EGR2"
## [1] "SNPs_NOTCH1"
## [1] "SNPs_POT1"
## [1] "SNPs_SF3B1"
## [1] "SNPs_TP53"
## [1] "SNPs_XPO1"
## [1] "health_record_bin_IGHV_mutated"
## [1] "health_record_bin_treated"limma_results$mut %>% unique## [1] "chrom_abber_del11q" "chrom_abber_del13q14"
## [3] "chrom_abber_del17p13" "chrom_abber_gain8q24"
## [5] "chrom_abber_trisomy12" "SNPs_ATM"
## [7] "SNPs_BIRC3" "SNPs_EGR2"
## [9] "SNPs_NOTCH1" "SNPs_POT1"
## [11] "SNPs_SF3B1" "SNPs_TP53"
## [13] "SNPs_XPO1" "health_record_bin_IGHV_mutated"
## [15] "health_record_bin_treated"Plot limma results
ggplot(data = limma_results) +
geom_histogram(aes(pvalue.limma, alpha = 0.5), bins = 40) +
guides(alpha = FALSE) +
xlab("p-value") +
facet_wrap( ~ mut +comparison, scale = "free_y") +
coord_cartesian(xlim = c(0, 1)) +
pp_sra + ggtitle("p-value histograms")
ggplot(data = limma_results) +
geom_histogram(aes(pvalue.limma_beforeDEQMS, alpha = 0.5), bins = 40) +
guides(alpha = FALSE) +
xlab("p-value") +
facet_wrap( ~ mut +comparison, scale = "free_y") +
coord_cartesian(xlim = c(0, 1)) +
pp_sra + ggtitle("p-values before DEqMS")
Limma Annotation of hits
Annotate
limma_results <- Annotate_Limma_Results(limma_results)
with(limma_results, table(mut, hit_annotation))## hit_annotation
## mut hit candidate no hit
## chrom_abber_del11q 2 6 7303
## chrom_abber_del13q14 0 3 7308
## chrom_abber_del17p13 9 14 7288
## chrom_abber_gain8q24 1 0 7310
## chrom_abber_trisomy12 205 332 6774
## health_record_bin_IGHV_mutated 48 211 7052
## health_record_bin_treated 3 72 7236
## SNPs_ATM 0 0 7311
## SNPs_BIRC3 0 4 7307
## SNPs_EGR2 4 10 7297
## SNPs_NOTCH1 0 0 7311
## SNPs_POT1 0 0 7311
## SNPs_SF3B1 61 43 7207
## SNPs_TP53 4 25 7282
## SNPs_XPO1 2 1 7308message("Number of up and downregulated hits")## Number of up and downregulated hitswith(limma_results %>% filter(hit==TRUE), table(mut, sign(logFC)))##
## mut -1 1
## chrom_abber_del11q 2 0
## chrom_abber_del17p13 8 1
## chrom_abber_gain8q24 1 0
## chrom_abber_trisomy12 48 157
## health_record_bin_IGHV_mutated 21 27
## health_record_bin_treated 3 0
## SNPs_EGR2 0 4
## SNPs_SF3B1 53 8
## SNPs_TP53 0 4
## SNPs_XPO1 2 0with(limma_results, table(mut, hit)) %>% as_tibble() %>%
filter(mut!="health_record_bin_treated") %>%
mutate(mut=gsub("chrom_abber", "CNVs", mut), mut=gsub("health_record_bin", "IGHV", mut)) %>%
separate(mut, into = c("Type", "mut"), sep="_", extra = "merge" ) %>%
filter(hit==TRUE) %>%
arrange(desc(n)) %>%
mutate(mut = as.factor(mut)) %>%
mutate(mut = factor(mut, levels = .$mut)) %>%
ggplot(aes(mut, n)) +
geom_col() +
ylab("Number of differentially abundant proteins") +
facet_grid(~Type, scales = "free_x", space="free") +
pp_sra_noguides_tilted
dif_proteins_P_plot <- with(limma_results %>% mutate(dir=sign(logFC)), table("mut"=paste0(mut,"dir", dir), hit)) %>% as_tibble() %>%
separate(mut, into = c("mut", "dir"), sep = "dir") %>%
filter(mut!="health_record_bin_treated") %>%
mutate(dir=as.numeric(dir)) %>%
mutate(n=n*dir) %>%
mutate(mut=gsub("chrom_abber", "CNVs", mut), mut=gsub("health_record_bin", "IGHV", mut)) %>%
separate(mut, into = c("Type", "mut"), sep="_", extra = "merge" ) %>%
filter(hit==TRUE) %>%
arrange(desc(n)) %>%
mutate(mut = as.factor(mut)) %>%
mutate(mut = factor(mut, levels = unique(.$mut))) %>%
ggplot(aes(mut, n)) +
geom_col(aes(fill=as.character(dir))) +
ylab("Nr. differentially abundant proteins") +
facet_grid(~Type, scales = "free_x", space="free") +
pp_sra_noguides_tilted +
scale_fill_manual(values = c("#0571b0", "#ca0020"), drop=FALSE) +
geom_hline(yintercept = 0, color="darkgray")
dif_proteins_P_plot
Plot limma hits in volcano plot
All conditions
ggplot(data = limma_results, aes(logFC, -log10(pvalue), colour = hit_annotation)) +
geom_vline(aes(xintercept = 0)) +
geom_point() +
geom_text(aes(label = rowname),
data = subset(limma_results, hit_annotation == "hit"),
vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
facet_wrap( ~ mut , ncol = 2) +
xlab("log2(fold change)") +
scale_colour_hue() +
pp_sra
tmp <- limma_results %>% filter(hit==TRUE) %>%
mutate(direction=sign(logFC)) %>%
group_by(mut) %>% arrange(fdr) %>%
dplyr::select(rowname, direction)## Adding missing grouping variables: `mut`message("Upregulated hits")## Upregulated hitssapply(group_nest(tmp)$mut, function(m){
tmp %>% filter(mut==m, direction==1) %>% .$rowname
})## $chrom_abber_del11q
## character(0)
##
## $chrom_abber_del17p13
## [1] "REXO4"
##
## $chrom_abber_gain8q24
## character(0)
##
## $chrom_abber_trisomy12
## [1] "SCYL2" "EEA1" "GOLGA3" "DERA" "TFCP2" "EIF4B"
## [7] "PTPN11" "CAMKK2" "CS" "TIGAR" "ANKLE2" "ESYT1"
## [13] "RAB5B" "ARID2" "EIF2B1" "ITPRIP" "GCN1" "TBC1D15"
## [19] "BRAP" "DENR" "USP15" "PSMD9" "PTPN6" "LDHB"
## [25] "IRAK4" "HIP1R" "TRIM25" "PYCARD" "OSBPL8" "TSFM"
## [31] "NAP1L1" "STRAP" "RAP1B" "PPP1R12A" "PYM1" "PBRM1"
## [37] "PA2G4" "NUDT1" "PRICKLE1" "PNP" "C12orf10" "MLX"
## [43] "PGAM5" "SAMHD1" "GNB4" "NR2C1" "SLC2A6" "RPAP3"
## [49] "CUL3" "SHMT2" "OGFOD2" "MON2" "RILPL2" "PIKFYVE"
## [55] "BTG1" "CORO1C" "RAN" "SART3" "TRAFD1" "PUS1"
## [61] "TWF1" "WNK1" "ARHGDIB" "CAND1" "RHOF" "XPOT"
## [67] "USP5" "MCTS1" "FRYL" "HECTD4" "TXNRD1" "CCDC91"
## [73] "STAT2" "RAB35" "ALDH16A1" "AGAP2" "TAOK3" "PARP12"
## [79] "APPL2" "IKBIP" "BRD7" "SPAG1" "SCAF11" "GRPEL1"
## [85] "ULK1" "CKLF" "PPP2R5B" "NT5E" "BCAT1" "CLIP2"
## [91] "ITGAL" "ACAD10" "FRMD8" "SAMD9" "ATXN2" "UBE2H"
## [97] "ABCD2" "MADD" "NLRC4" "IPO13" "SCPEP1" "MLXIP"
## [103] "SLC25A3" "ASPSCR1" "BCL7A" "ERP29" "TTC39C" "AKT2"
## [109] "ATP6V0A2" "METAP2" "FES" "CCDC92" "NCOR2" "SETD1B"
## [115] "RAB21" "CYFIP1" "FLNA" "CNOT8" "RASSF3" "LDB1"
## [121] "TMEM19" "CALCOCO1" "ANKRD13A" "ITGB2" "DPYSL2" "DNAJC7"
## [127] "ZDHHC17" "PIP4K2C" "PRKAB1" "TMEM231" "KNTC1" "EFHD2"
## [133] "RAB29" "RHBDF2" "INF2" "FIBP" "PIH1D1" "C2CD5"
## [139] "STAT6" "METTL25" "RNF135" "MAP3K8" "NCKAP5L" "BIN2"
## [145] "ITGAX" "SZRD1" "TMEM102" "STK38" "WASHC3" "ATF7"
## [151] "PITPNA" "UNG" "HCFC2" "EEF1A1" "SH3GL1" "YARS2"
## [157] "ATP2A2"
##
## $health_record_bin_IGHV_mutated
## [1] "FTH1" "FTL" "USP6NL" "GRK3" "TBC1D2B" "ZBTB20" "ABI3"
## [8] "FCRL1" "SLAMF1" "FCRL2" "UBE2E2" "ADD3" "DUSP22" "ADD1"
## [15] "MTSS1" "SMC6" "GNA13" "KLHL14" "ELMO2" "PLEKHF2" "SMC5"
## [22] "FCHSD2" "PIP4K2A" "RASA2" "FLYWCH2" "STK17B" "POT1"
##
## $health_record_bin_treated
## character(0)
##
## $SNPs_EGR2
## [1] "SLC2A1" "ACTN1" "EPPK1" "ZNF446"
##
## $SNPs_SF3B1
## [1] "SUGP1" "PSPC1" "SMARCAD1" "DGCR8" "ATAD3B" "EYA3" "SLC4A1AP"
## [8] "MSH2"
##
## $SNPs_TP53
## [1] "LRPAP1" "TP53" "REXO4" "FUT11"
##
## $SNPs_XPO1
## character(0)message("Downregulated hits")## Downregulated hitssapply(group_nest(tmp)$mut, function(m){
tmp %>% filter(mut==m, direction==-1) %>% .$rowname
})## $chrom_abber_del11q
## [1] "CUL5" "AASDHPPT"
##
## $chrom_abber_del17p13
## [1] "FXR2" "MAP2K4" "PAFAH1B1" "GLOD4" "ZBTB4" "ELAC2" "YWHAE"
## [8] "RABEP1"
##
## $chrom_abber_gain8q24
## [1] "MCPH1"
##
## $chrom_abber_trisomy12
## [1] "FRA10AC1" "PELI1" "TGFBR2" "TCEAL3" "BCL7C" "ABHD14B"
## [7] "GDPD1" "NUMA1" "SRSF5" "RANBP3" "ACYP1" "ZFP64"
## [13] "PARP2" "CRAT" "TCF4" "TLK1" "BCOR" "WDR89"
## [19] "PTPRA" "SELENOM" "TLE4" "DNTTIP1" "CAMK2D" "MORF4L2"
## [25] "PCBD1" "ELMSAN1" "NOL7" "GPATCH8" "UBXN7" "NCOR1"
## [31] "PLEKHG1" "PPP4R3A" "IRS1" "SAMSN1" "TCEAL4" "ROR1"
## [37] "C3orf38" "LANCL1" "CEP44" "GPD2" "ANP32E" "CTH"
## [43] "MBD2" "TTC33" "DBR1" "TFAP4" "SLC25A23" "VAT1"
##
## $health_record_bin_IGHV_mutated
## [1] "CPT1A" "TMEM165" "GALNT2" "UBE2H" "FAM114A1" "SERPINB6"
## [7] "UBAP2" "ATP2C1" "VAMP5" "SOGA1" "ZAP70" "CCDC186"
## [13] "ATOX1" "ATP2A2" "SLC35B2" "APOL3" "PDXDC1" "PLCH1"
## [19] "MFN1" "SEPT10" "APOBEC3G"
##
## $health_record_bin_treated
## [1] "MTSS1" "ICAM2" "TNIP1"
##
## $SNPs_EGR2
## character(0)
##
## $SNPs_SF3B1
## [1] "TPP2" "MAP3K7" "PPP2R5A" "WASHC5" "WASHC4" "VPS8"
## [7] "IQGAP1" "WASHC1" "NAA16" "FKBP15" "DPH5" "TAB1"
## [13] "TTI1" "CCDC88B" "SEPSECS" "TRAPPC6B" "SEPT6" "TTI2"
## [19] "WASHC3" "ABCB7" "DDX58" "INPPL1" "EML3" "SAFB2"
## [25] "YWHAB" "WASHC2A" "C12orf66" "TTC37" "SKIV2L" "DIP2A"
## [31] "PDS5A" "NPRL2" "IDUA" "RECQL" "UQCC1" "TNPO3"
## [37] "CEP135" "SZT2" "PPP6R3" "AP2A2" "ARIH1" "WASHC2C"
## [43] "JMY" "C10orf76" "LIG3" "SUN1" "SCLT1" "DEPDC5"
## [49] "MTMR6" "KPNA1" "PNKP" "MROH1" "SEPT2"
##
## $SNPs_TP53
## character(0)
##
## $SNPs_XPO1
## [1] "XPO1" "DHX36"Trisomy 12 colored by chromosome
limma_results %>% filter(mut=="chrom_abber_trisomy12") %>%
mutate("Affected_region" = if_else( chromosome_name=="12", "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
mutate(logFC= if_else(logFC > log2(5), log2(5.2),
if_else(logFC<(-log2(5)), -log2(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname, color=chromosome_name=="12"),
data = subset(limma_results %>% filter(mut=="chrom_abber_trisomy12"), hit_annotation == "hit"),
vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
guides(color=FALSE)+
coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue)))) +
ggtitle("trisomy12")+
theme( plot.title = element_text(size = 20)) 
tmp <- limma_results %>% filter(mut=="chrom_abber_trisomy12", chromosome_name=="12") %>%
.$hit_annotation %>% table %>% as_tibble()
colnames(tmp) <- c("Protein", "n")
tmptmp %>%
mutate(Protein=as.factor(Protein)) %>%
mutate(Protein=factor(Protein, levels = c("hit", "candidate", "no hit"))) %>%
ggplot(aes(x="", fill=Protein,y= n)) +
geom_col() +
ggtitle("Hit annotation proteins on chromosome 12 in trisomy 12") +
coord_polar("y", start=0) +
pp_sra +
theme(axis.title = element_blank(), axis.ticks = element_blank())
tmp <- (limma_results %>%
filter(mut=="chrom_abber_trisomy12", chromosome_name=="12") %>%
.$logFC>0) %>% table %>% as_tibble()
colnames(tmp) <- c("fc", "n")
tmp <- tmp %>% arrange(desc(`fc`))
tmp <- tmp %>% mutate(label=paste( round(n/sum(tmp$n)*100, 1), "%" ), cums =cumsum(tmp$n) )
tmp %>%
mutate("fc"=as.factor(`fc`)) %>%
ggplot(aes(x="", fill=`fc`,y= n)) +
geom_col() +
ggtitle("Proteins on chromosome 12 with positive logFC in trisomy 12") +
coord_polar("y", start=0) +
pp_sra +
theme(axis.title = element_blank(), axis.ticks = element_blank())+
scale_fill_manual(values = c( "grey", "#0571b0"), drop=FALSE)+
annotate(geom = "text", y = tmp$cums-(tmp$n/2), x = 1, label = tmp$label)
tmp <- (limma_results %>%
filter(mut=="chrom_abber_trisomy12", hit==TRUE) %>%
.$chromosome_name == "12") %>%
table %>% as_tibble()
colnames(tmp) <- c("chr12", "n")
tmp <- tmp %>% arrange(desc(`chr12`))
tmp <- tmp %>% mutate(label=paste( round(n/sum(tmp$n)*100, 1), "%" ), cums =cumsum(tmp$n) )
piechart_hits_prot_chr_tris12 <- tmp %>%
mutate("chr12"=as.factor(`chr12`)) %>%
ggplot(aes(x="", fill=`chr12`,y= n)) +
geom_col() +
ggtitle("Chromosome location of hits in trisomy 12") +
coord_polar("y", start=0) +
pp_sra +
theme(axis.title = element_blank(), axis.ticks = element_blank())+
scale_fill_manual(values = c( "grey", "#0571b0"), drop=FALSE)
piechart_hits_prot_chr_tris12 +
annotate(geom = "text", y = tmp$cums-(tmp$n/2), x = 1, label = tmp$label)
del13q14 colored by chromosome
limma_results %>% filter(mut=="chrom_abber_del13q14") %>%
mutate("Affected_region" = if_else( (chromosome_name=="13" & start_position>47000000 & start_position<51000000 ), "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
mutate(logFC= if_else(logFC > log2(5), log2(5.2),
if_else(logFC<(-log2(5)), -log(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname, color=(chromosome_name=="13" & start_position>47000000 & start_position<51000000)),
data = subset(limma_results %>% filter(mut=="chrom_abber_del13q14"), hit_annotation == "hit"),
vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
guides(color=FALSE)+
coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue))))+
ggtitle("del13q14")+
theme( plot.title = element_text(size = 20)) 
del11q colored by chromosome
limma_results %>% filter(mut=="chrom_abber_del11q") %>%
mutate("Affected_region" = if_else( (chromosome_name=="11" & start_position>97400000 & start_position<110600000), "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
mutate(logFC= if_else(logFC > log2(5), log2(5.2),
if_else(logFC<(-log2(5)), -log2(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname, color=(chromosome_name=="11" & start_position>97400000 & start_position<110600000)),
data = subset(limma_results %>% filter(mut=="chrom_abber_del11q"), hit_annotation == "hit"),
vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
guides(color=FALSE)+
coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue)))) +
ggtitle("del11q")+
theme(plot.title = element_text(size = 20)) 
del17p13 colored by chromosome
limma_results %>% filter(mut=="chrom_abber_del17p13") %>%
mutate("Affected_region" = if_else( (chromosome_name=="17" & mean_position<10800000), "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
mutate(logFC= if_else(logFC > log2(5), log2(5.2),
if_else(logFC<(-log2(5)), -log2(5.2), logFC) ) ) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname, color=(chromosome_name=="17" & mean_position<10800000)),
data = subset(limma_results %>% filter(mut=="chrom_abber_del17p13"), hit_annotation == "hit"),
vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
guides(color=FALSE)+
coord_cartesian(xlim=c(-log2(5), log2(5)), ylim=c(0, -log10(min(limma_results$pvalue)))) +
ggtitle("del17p13")+
theme(plot.title = element_text(size = 20)) 
XPO1
XPO_volcano_plot <- limma_results %>% filter(mut=="SNPs_XPO1") %>%
mutate("Affected_region" = if_else( rowname=="XPO1", "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname),
data = subset(limma_results %>% filter(mut=="SNPs_XPO1", rowname!="XPO1"), hit_annotation == "hit"),
vjust = 0, size = 3, nudge_y = (-1), nudge_x = -0.3,# nudge_y = -0.3,
check_overlap = FALSE, color="#0571b0") +
geom_text(aes(label = rowname, color=rowname=="XPO1"),
data = subset(limma_results %>% filter(mut=="SNPs_XPO1", rowname=="XPO1")),
vjust = 0, size = 3, nudge_x = 0.2, nudge_y = (-0.5), #nudge_y = (-0.1),
check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
coord_cartesian(xlim=c(-log2(5), log2(5)) )
XPO_volcano_plot + ggtitle("XPO1")+
theme( plot.title = element_text(size = 20)) +
guides(color=FALSE)
TP53
TP53_volcano_plot <- limma_results %>% filter(mut=="SNPs_TP53") %>%
mutate("Affected_region" = if_else( rowname=="TP53", "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname),
data = subset(limma_results %>% filter(mut=="SNPs_TP53", rowname!="TP53"), hit_annotation == "hit"),
vjust = 0, size = 3, nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE, color="#0571b0") +
geom_text(aes(label = rowname, color=rowname=="TP53"),
data = subset(limma_results %>% filter(mut=="SNPs_TP53", rowname=="TP53")),
vjust = 0, size = 3, nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
coord_cartesian(xlim=c(-log2(2.2), log2(2.2)) )
TP53_volcano_plot + ggtitle("TP53")+
theme( plot.title = element_text(size = 20)) +
guides(color=FALSE)
ATM
ATM_volcano_plot <- limma_results %>% filter(mut=="SNPs_ATM") %>%
mutate("Affected_region" = if_else( rowname=="ATM", "TRUE",
if_else(hit_annotation == "hit", "FALSE", "not altered"))) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("FALSE","TRUE", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname),
data = subset(limma_results %>% filter(mut=="SNPs_ATM", rowname!="ATM"), hit_annotation == "hit"),
vjust = 0, size = 3, nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE, color="#0571b0") +
geom_text(aes(label = rowname, color=rowname=="ATM"),
data = subset(limma_results %>% filter(mut=="SNPs_ATM", rowname=="ATM")),
vjust = 0, size = 3, nudge_y = -0.3, nudge_x = -0.2, check_overlap = FALSE) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "orange1", "grey"), drop=FALSE) +
pp_sra+
coord_cartesian(xlim=c(-log2(2.2), log2(2.2)) )
ATM_volcano_plot + ggtitle("ATM")+
theme( plot.title = element_text(size = 20)) +
guides(color=FALSE)
Limma with subsets of patients
Consensus clusters groups 5 vs. all others
selpats <- colData(multiomics_MAE) %>% as_tibble() %>% filter(!is.na(PG)) %>% .$patient_ID
proteomics_tbl_meta_biomart_tmp <- proteomics_tbl_meta_biomart
proteomics_tbl_meta_biomart <- left_join(proteomics_tbl_meta_biomart, enframe(CCP_group5==5, name = "primary", value="PG"), by="primary" )
limma_results_CC5_5_all <- Calculate_Limma(assay_data = assay(multiomics_MAE[prot_few_nas , selpats ,"proteomics"]), mut = "PG")## ExperimentList contains data.frame or DataFrame,
## potential for errors with mixed data types
## ExperimentList contains data.frame or DataFrame,
## potential for errors with mixed data typesproteomics_tbl_meta_biomart <- proteomics_tbl_meta_biomart_tmp
ggplot(data = limma_results_CC5_5_all) +
geom_histogram(aes(pvalue.limma, alpha = 0.5), bins = 40) +
guides(alpha = FALSE) +
xlab("p-value") +
facet_wrap( ~ mut , scale = "free_y") +
coord_cartesian(xlim = c(0, 1)) +
pp_sra +
ggtitle("PG groups")
limma_results_CC5_5_all <- Annotate_Limma_Results(limma_results_CC5_5_all)
limma_results_CC5_5_all %>%
mutate("Affected_region" = if_else(hit_annotation == "hit", "hit", "not altered")) %>%
mutate("Affected_region"= factor(Affected_region, levels=c("hit", "not altered"))) %>%
ggplot(aes(logFC, -log10(pvalue))) +
geom_vline(aes(xintercept = 0)) +
geom_point(alpha=0.8, aes(colour = Affected_region)) +
geom_text(aes(label = rowname), color="#0571b0",
data = subset(limma_results_CC5_5_all %>% filter( hit_annotation == "hit"),
vjust = 0, nudge_y = 0.1, size = 3, check_overlap = FALSE)) +
facet_wrap( ~ mut, ncol = 1) +
xlab("log2(fold change)") +
scale_color_manual(values = c("#0571b0", "grey"), drop=FALSE) +
pp_sra+
guides(color=guide_legend(title = "PG groups"))+
ggtitle("Consensus cluster plus group 5 against all others") 
limma_results_CC5_5_all %>% filter(hit==TRUE) %>% .$rowname## [1] "WDR1" "ACTR2" "ACTR3" "FNTA" "CAP1" "ATG5"
## [7] "SH3BGRL" "GAPVD1" "NRBP1" "CDC5L" "PITPNB" "BTK"
## [13] "ARPC3" "TERF2" "ARPC1B" "SARS" "LUC7L3" "ATG7"
## [19] "VAV2" "CFL1" "ARPC5" "IKBKB" "NUP155" "PGK1"
## [25] "TBC1D7" "CWC25" "DHX9" "PPP6C" "PCF11" "PANK4"
## [31] "PIK3CD" "NUP98" "TERF2IP" "AFP" "NUP54" "SRSF4"
## [37] "ARPC2" "RALGAPB" "PDP1" "NUP93" "PUF60" "OTULIN"
## [43] "GLE1" "CSK" "GRK2" "NUP133" "RANBP1" "ARHGDIA"
## [49] "USP9X" "ARHGAP17" "ANKRD44" "SNX5" "SF3B1" "RABGGTB"
## [55] "UBR1" "C16orf58" "HNRNPC" "ZNRF2" "HTT" "ARFGEF1"
## [61] "SF3B6" "PNN" "CRNKL1" "SFPQ" "GSR" "CORO1A"
## [67] "NUP58" "PIK3R1" "RBM12B" "METAP2" "TWF2" "SH3GL1"
## [73] "OSTF1" "NUP205" "MOB2" "YWHAB" "TLN1" "PPP6R1"
## [79] "SLTM" "OTUB1" "EMC3" "PAXBP1" "HNRNPR" "SMPD4"
## [85] "STAT6" "MTX1" "PGGT1B" "PRPF19" "DDX46" "WIPF1"
## [91] "TRAPPC10" "USP47" "PRKAG1" "CCDC47" "NUP210" "MAPRE1"
## [97] "EEF2" "ASXL2" "RPN2" "BUD13" "RAC2" "UBE2O"
## [103] "CRKL" "CAPZB" "GDI2" "MTMR9" "RPN1" "SNRNP48"
## [109] "PAF1" "KDELC1" "NXF1" "THRAP3" "HECTD1" "SON"
## [115] "CIAO1" "EEF1A1" "TBCA" "RASSF2" "AKT2" "ATG3"
## [121] "SAP18" "MAP2K2" "SERPINC1" "UBE2V2" "DENND1C" "YEATS2"
## [127] "NONO" "EFHD2" "FBXW2" "NDUFAF2" "PAPOLG" "YLPM1"
## [133] "SUPT6H" "SSR4" "MARK2" "UBE2N" "TNFAIP8" "BIN2"
## [139] "HAX1" "PTGES3" "SRSF3" "DIAPH1" "SF3B2" "IMMT"
## [145] "IST1" "UBAC1" "ZFR" "DCTN4" "MTMR6" "YKT6"
## [151] "ZNF512" "CPSF3" "KHDRBS1" "PTPA" "GTPBP1" "CDC73"
## [157] "PDPR" "TF" "TSC1" "CCDC25" "WASF2" "BRK1"
## [163] "HMGXB4" "GRB2" "RSRC1" "CLEC3B" "MYD88" "PLEKHM2"
## [169] "INPP5D" "NUP107" "HNRNPM" "ZBTB48" "C16orf70" "GFM1"
## [175] "TPD52L2" "RBBP6" "PTPN11" "SNRNP70" "FBLN1" "PDCD6IP"
## [181] "CAPZA1" "YTHDC1" "NFKBIE" "DARS2" "CYC1" "DDHD2"
## [187] "COMMD10" "PFN1" "NUP62" "MYL6" "PSMB2" "NCKAP1L"
## [193] "EMC1" "STAT5A" "HM13" "RNPS1" "UVRAG" "PACS1"
## [199] "ILF3" "XPNPEP1" "PIK3C3" "INPPL1" "CARS" "PNPLA8"
## [205] "NFKB1" "DDX10" "DBNL" "NAP1L4" "CERS2" "TOMM70"
## [211] "SF3A3" "VPS4B" "SF3A1" "RBM17" "BAP1" "GNPDA2"
## [217] "TRAPPC1" "EXOC4" "ARSA" "CAMLG" "GSPT1" "POLE4"
## [223] "GLS" "MBIP" "NUP88" "NUPL2" "LRBA" "DCTN6"
## [229] "APMAP" "CHMP4B" "UTP3" "MAP2K1" "POMGNT1" "GMFB"
## [235] "SPCS2" "LTA4H" "RDH14" "PRKAA1" "RALGAPA2" "GAK"
## [241] "CENPB" "KYAT1" "RBM26" "CYFIP2" "SLC25A12" "UBA5"
## [247] "NRAS" "USP24" "DDX50" "PPP3R1" "IFI16" "ADAR"
## [253] "REL" "ATP6V1C1" "NUP188" "PCCB" "MSL1" "SRRM2"
## [259] "RMC1" "RABL6" "RIF1" "ELP2" "RNF6" "EMC2"
## [265] "WDR37" "SUZ12" "TMEM214" "RIC1" "SRSF11" "HIBADH"
## [271] "TRAPPC9" "PREP" "PDPK1" "RANBP2" "PPP3CB" "VAPA"
## [277] "CDK14" "DOCK11" "TSC2" "PRPSAP1" "OFD1" "DKC1"
## [283] "WAS" "PIBF1" "OCIAD1" "EXOC2" "RABGGTA" "ALDH5A1"
## [289] "CTR9" "TPR" "HNRNPH2" "HSH2D" "METTL21A" "FAM160B1"
## [295] "GOLGA5" "RNF126" "DOCK2" "ESF1" "PLA2G15" "DCTN2"
## [301] "RBM39" "SSR3" "ACAP2" "RBM15" "ATP6V1B2" "IARS2"
## [307] "ITIH2" "CCT7" "BRAP" "RPS12" "DNMT1" "TCP1"
## [313] "PLCG2" "MOSPD2" "AKAP13" "ATP6V1H" "DAPP1" "TFIP11"
## [319] "TMX1" "ZFC3H1" "CLPB" "SNRPA1" "CISD2" "CNN2"
## [325] "MALT1" "MAPK1" "HABP2" "ENOPH1" "DOK1" "PCCA"
## [331] "MFN2" "PTP4A2" "SLU7" "CCAR1" "ANKRD13A" "PSMA6"
## [337] "EGLN1" "FLII" "COMMD9" "ELMOD2" "SAFB" "PKN1"
## [343] "FKBP1A" "CYB5B" "PRPS1" "CAND1" "NOL9" "PGK2"
## [349] "TRMT61B" "NTMT1" "GGCT" "VAV1" "SNRPB2" "ATP7A"
## [355] "HCLS1" "NELFCD" "HACD3" "PDLIM2" "ARAF" "PTPN6"
## [361] "PDCD7" "CLEC16A" "MLYCD" "PDK3" "PES1" "GPN1"
## [367] "ARMT1" "SNW1" "CBX3" "SRSF5" "AKT1" "BECN1"
## [373] "SRSF10" "MAPK14" "CENPC" "CPSF2" "PBDC1" "TECPR1"
## [379] "APPL1" "A2M" "SLC25A36" "PLAA" "AIFM1" "ZMPSTE24"
## [385] "ROCK1" "MAPK15" "GMFG" "TMEM208" "BAZ1B" "EMC7"
## [391] "ETF1" "PTEN" "TRAM1" "PDHA1" "HOOK3" "RABGAP1"
## [397] "FUNDC2" "VAMP4" "MCM3AP" "SMARCA5" "TRADD" "MSN"
## [403] "TST" "IKBKG" "RMDN1" "UNC45A" "RFX1" "APAF1"
## [409] "ATG16L1" "IK" "ATP6V1A" "WDR44" "CLASRP" "NUP35"
## [415] "RPL12" "CSTF3" "TRIP11" "ABL1" "CANX" "NIPSNAP2"
## [421] "STIM2" "UBR2" "PAWR" "ISYNA1" "PLEKHF2" "PLS1"
## [427] "SEC22B" "KMT2B" "TMCC1" "RABAC1" "DDRGK1" "KIAA0100"
## [433] "CAPN1" "ABI1" "TTC9" "GNL1" "TFAM" "STT3B"
## [439] "GLTP" "VPS33A" "RASGRP2" "HMGN3" "MTMR12" "IRAK4"
## [445] "ZC3H18" "PDHB" "PAAF1" "DCTN5" "GAR1" "SPTY2D1"
## [451] "DHX16" "VPS37A" "GOT1" "WTAP" "LIMD1" "WDR5"
## [457] "AP3M1" "ELP4" "CLMN" "SPCS3" "RSBN1L" "RPL14"
## [463] "RHOT1" "SLC25A32" "ITCH" "DDOST" "IKZF1" "DNAJC11"
## [469] "CHCHD6" "SFXN3" "TNIP2" "CLIC1" "EML3" "ARHGDIB"
## [475] "ERGIC1" "EXOC3" "CEP44" "AUP1" "EXOC1" "RNF123"
## [481] "MBLAC2" "CAPNS1" "RPE" "ALDH18A1" "HTATIP2" "RAD23A"
## [487] "VDAC1" "VCPIP1" "BPTF" "DENND4C" "GPATCH4" "SF3B4"
## [493] "CCDC90B" "UBAP1" "IWS1" "RBM10" "CTNNBL1" "TARS"
## [499] "NDC1" "PIK3R4" "DRG2" "DIDO1" "CCT8" "SLC38A2"
## [505] "U2SURP" "BCL10" "RIC8A" "SF1" "MTPN" "BLOC1S3"
## [511] "MAPKAPK5" "GATAD1" "MED1" "SRSF7" "PTK2B" "MUL1"
## [517] "EIF5" "GSTCD" "MTMR14" "PLRG1" "NECAP2" "STAT3"
## [523] "XRN2" "BLOC1S5" "PIGX" "SSRP1" "CHCHD3" "EMC4"
## [529] "CHD4" "LSG1" "RACK1" "FBXO21" "SYMPK" "PPIA"
## [535] "CUL3" "ZNF638" "FAM76A" "SH2D3C" "COTL1" "TBCB"
## [541] "HADHB" "MORC2" "BRAF" "CHMP2B" "ATP6AP2" "SMC3"
## [547] "XAB2" "CLPTM1L" "MPHOSPH10" "NELFA" "KYNU" "SRSF1"
## [553] "EXOC7" "LMAN2" "MTIF3" "NFATC2" "PSMB9" "RAD21"
## [559] "AKT3" "SEC61B" "MCU" "MTDH" "SRRM1" "YWHAG"
## [565] "FGD3" "ELMO2" "BLOC1S2" "C2CD2L" "TRAPPC6B" "OGFR"
## [571] "PDZD8" "MAD1L1" "ILF2" "RALGPS2" "ELP3" "TXNRD1"
## [577] "EPM2A" "UBE2Z" "MYL12A" "OLA1" "FUNDC1" "RETSAT"
## [583] "PREB" "CDC42EP3" "STXBP2" "RPRD2" "BLOC1S4" "DCTN3"
## [589] "LRRFIP1" "CDK7" "SNRPD2" "LCP1" "CLPX" "ATG4B"
## [595] "CDC42" "AKAP1" "GEMIN5" "GFM2" "KIAA0753" "PFKL"
## [601] "PI4KB" "CAB39" "CCT3" "NKIRAS2" "DCTN1" "GON4L"
## [607] "UBE2M" "RABGEF1" "TRMT1L" "DUT" "MRPL15" "NAB1"
## [613] "UBA6" "RGL2" "GPATCH8" "ANKMY2" "DTNBP1" "PPP1R12C"
## [619] "CAPN7" "RPL30" "VMP1" "COPS7B" "GATAD2B" "TRAPPC3"
## [625] "ZC3H4" "FERMT3" "SIN3B" "PPP3CC" "RPSA" "KPNA2"
## [631] "ZCCHC6" "MYO1E" "GOLIM4" "ATAD2" "VPS35L" "HADH"
## [637] "KIF3B" "SPOP" "BIN1" "PRKAB1" "PRKD2" "RPL22"
## [643] "TMED10" "PSMD8" "MRGBP" "PRKCZ" "RRAS2" "ATP6V1G1"
## [649] "MSANTD4" "AHCTF1" "YTHDF3" "ARPC5L" "UBR4" "PTPN7"
## [655] "CEP162" "SNX2" "CBFB" "NARS" "STT3A" "NSDHL"
## [661] "CDKN2AIP" "KIF5B" "ACOX1" "SLC22A18" "ARIH1" "NIPSNAP1"
## [667] "MCUB" "NAP1L1" "RP2" "KIFAP3" "SREK1" "TPP1"
## [673] "UFL1" "PPOX" "PPP1R16B" "STAP1" "ELP1" "CNOT9"
## [679] "NDUFAF6" "TCERG1" "VIRMA" "TARS2" "MDH2" "GRK3"
## [685] "NUP160" "SERPINF2" "USP25" "CLIC4" "RBMX" "RRP9"
## [691] "RRP1" "PRDX3" "IMP3" "AIMP2" "BTF3L4" "PHPT1"
## [697] "TBC1D9B" "CNDP2" "SSBP1" "WASHC3" "PRAG1" "FMC1"
## [703] "TTI1" "PDHX" "TOM1" "NSMAF" "NAPEPLD" "GDI1"
## [709] "MAPK9" "ZNF326" "MICU1" "RIPOR2" "NCOA5" "NLN"
## [715] "SNX3" "LYAR" "ZBTB9" "NUP50" "FNBP4" "C3orf33"
## [721] "GPN2" "DOPEY2" "PSMA5" "ADRM1" "FAF2" "SH3BP2"
## [727] "ZNHIT1" "MIEN1" "CHRAC1" "SMARCE1" "CCDC43" "PIGK"
## [733] "ERCC4" "MOB3A" "KNOP1" "RPL38" "EXOC8" "PXK"
## [739] "XXYLT1" "SLC4A1AP" "HMGCL" "MTX2" "LIMD2" "TMED2"
## [745] "MAN1B1" "MYO9B" "HNRNPH3" "MAP3K8" "USP3" "RIPK1"
## [751] "IGBP1" "GNB2" "ATP6V1E1" "CIRBP"Gene set enrichment analysis consensus cluster plus groups 5 vs. all
####### Input of all Values
longFormat((multiomics_MAE[prot_few_nas , selpats ,"proteomics"])) %>%
as_tibble() %>%
dplyr::select("NAME"=rowname, "DESCRIPTION"=assay , primary, value) %>%
unique() %>%
spread(primary, value) #%>%
#write.table(file = "/Users/sophierabe/Desktop/PhD/Labor/Proteomics/CLL/GSEA_CLL_Proteomics/190927_GSEA_data_long_gradient_CCP_5vsall.txt", quote = FALSE, row.names = FALSE, sep = "\t", na="")
tmp_pats <- longFormat((multiomics_MAE[prot_few_nas , selpats ,"proteomics"])) %>%
as_tibble() %>%
dplyr::select("NAME"=rowname, primary, value) %>%
unique() %>%
spread(primary, value) %>%
dplyr::select(-1) %>%
colnames() %>%
enframe()
tmp_pats <- tmp_pats %>%
left_join(.,
enframe(multiomics_MAE$PG, value="PG", name="primary" ) %>% mutate("PG"=PG==5) ,
by=c("value"="primary") )
tmp_pats %>%
dplyr::select(PG) #%>% t() %>%
#write.table(file = "/Users/sophierabe/Desktop/PhD/Labor/Proteomics/CLL/GSEA_CLL_Proteomics/190927_GSEA_phenotypelabels_long_gradient__CCP_5vsall.cls", quote = FALSE, row.names = FALSE, col.names = FALSE, sep = " ", na="")
dim(tmp_pats)
sum(tmp_pats$PG==TRUE)Save data needed for further analysis
save(limma_results,
file = "/Volumes/sd17b003/Sophie/Analysis/CLL_Proteomics/CLL_Proteomics_final/Proteomics_Git/Robjects/CLL_Proteomics_LimmaProteomics.RData")Session Info
sessionInfo()## R version 3.5.2 (2018-12-20)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Mojave 10.14.6
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] DEqMS_1.0.1 cowplot_1.0.0
## [3] survminer_0.4.6 ConsensusClusterPlus_1.46.0
## [5] RColorBrewer_1.1-2 gplots_3.0.1.1
## [7] PMA_1.1 MultiAssayExperiment_1.8.3
## [9] biomaRt_2.38.0 biomartr_0.9.0
## [11] Rtsne_0.15 IHW_1.10.1
## [13] apeglm_1.4.2 limma_3.38.3
## [15] fdrtool_1.2.15 vsn_3.50.0
## [17] forcats_0.4.0 stringr_1.4.0
## [19] dplyr_0.8.3 purrr_0.3.3
## [21] readr_1.3.1 tidyr_1.0.0
## [23] tibble_2.1.3 tidyverse_1.3.0
## [25] TOSTER_0.3.4 DESeq2_1.22.2
## [27] SummarizedExperiment_1.12.0 DelayedArray_0.8.0
## [29] BiocParallel_1.16.6 matrixStats_0.55.0
## [31] Biobase_2.42.0 GenomicRanges_1.34.0
## [33] GenomeInfoDb_1.18.2 IRanges_2.16.0
## [35] S4Vectors_0.20.1 BiocGenerics_0.28.0
## [37] readxl_1.3.1 ggrepel_0.8.1
## [39] ggpmisc_0.3.3 ggfortify_0.4.8
## [41] ggbeeswarm_0.6.0 effsize_0.7.6
## [43] ggpubr_0.2.4 magrittr_1.5
## [45] lemon_0.4.3 Hmisc_4.3-0
## [47] ggplot2_3.2.1 Formula_1.2-3
## [49] survival_3.1-8 lattice_0.20-38
## [51] Matrix_1.2-18 progress_1.2.2
## [53] reshape2_1.4.3 pheatmap_1.0.12
## [55] openxlsx_4.1.4 gtools_3.8.1
## [57] plyr_1.8.4
##
## loaded via a namespace (and not attached):
## [1] tidyselect_0.2.5 RSQLite_2.1.4 AnnotationDbi_1.44.0
## [4] htmlwidgets_1.5.1 munsell_0.5.0 codetools_0.2-16
## [7] preprocessCore_1.44.0 withr_2.1.2 colorspace_1.4-1
## [10] knitr_1.26 rstudioapi_0.10 ggsignif_0.6.0
## [13] labeling_0.3 slam_0.1-46 bbmle_1.0.20
## [16] GenomeInfoDbData_1.2.0 lpsymphony_1.10.0 KMsurv_0.1-5
## [19] farver_2.0.1 bit64_0.9-7 coda_0.19-3
## [22] vctrs_0.2.0 generics_0.0.2 xfun_0.11
## [25] R6_2.4.1 locfit_1.5-9.1 bitops_1.0-6
## [28] assertthat_0.2.1 scales_1.1.0 nnet_7.3-12
## [31] beeswarm_0.2.3 gtable_0.3.0 affy_1.60.0
## [34] rlang_0.4.2 zeallot_0.1.0 genefilter_1.64.0
## [37] splines_3.5.2 lazyeval_0.2.2 acepack_1.4.1
## [40] impute_1.56.0 broom_0.5.2 checkmate_1.9.4
## [43] BiocManager_1.30.10 yaml_2.2.0 modelr_0.1.5
## [46] backports_1.1.5 tools_3.5.2 ellipsis_0.3.0
## [49] affyio_1.52.0 Rcpp_1.0.3 base64enc_0.1-3
## [52] zlibbioc_1.28.0 RCurl_1.95-4.12 prettyunits_1.0.2
## [55] rpart_4.1-15 zoo_1.8-6 haven_2.2.0
## [58] cluster_2.1.0 fs_1.3.1 data.table_1.12.8
## [61] reprex_0.3.0 hms_0.5.2 evaluate_0.14
## [64] xtable_1.8-4 XML_3.98-1.20 emdbook_1.3.11
## [67] gridExtra_2.3 compiler_3.5.2 KernSmooth_2.23-16
## [70] crayon_1.3.4 htmltools_0.4.0 geneplotter_1.60.0
## [73] lubridate_1.7.4 DBI_1.0.0 dbplyr_1.4.2
## [76] MASS_7.3-51.4 cli_2.0.0 gdata_2.18.0
## [79] pkgconfig_2.0.3 km.ci_0.5-2 numDeriv_2016.8-1.1
## [82] foreign_0.8-72 xml2_1.2.2 annotate_1.60.1
## [85] vipor_0.4.5 XVector_0.22.0 rvest_0.3.5
## [88] digest_0.6.23 Biostrings_2.50.2 rmarkdown_1.18
## [91] cellranger_1.1.0 survMisc_0.5.5 htmlTable_1.13.3
## [94] curl_4.3 lifecycle_0.1.0 nlme_3.1-142
## [97] jsonlite_1.6 fansi_0.4.0 pillar_1.4.2
## [100] httr_1.4.1 glue_1.3.1 zip_2.0.4
## [103] bit_1.1-14 stringi_1.4.3 blob_1.2.0
## [106] latticeExtra_0.6-28 caTools_1.17.1.3 memoise_1.1.0